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The hydro-dynamic phase field model is applied to the 
problem of film spreading on a solid surface. The disjoin- 
ing potential, responsible for modification of the fluid prop- 
erties near a three-phase contact line, is computed from the 
solvability conditions of the density field equation with ap- 
propriate boundary conditions imposed on the solid support. 
The equation describing the motion of a spreading film are 
derived in the lubrication approximation. In the case of 
quasi-equilibrium spreading, is shown that the correct sharp- 
interface limit is obtained, and sample solutions are obtained 
by numerical integration. It is further shown that evapora- 
tion or condensation may strongly affect the dynamics near 
the contact line, and accounting for kinetic retardation of the 
interphase transport is necessary to build up a consistent the- 
ory. 



I. INTRODUCTION 

One of long-standing hydrodynamic riddles is the na- 
ture of viscous flow in the vicinity of a three-phase 
(gas-liquid-solid) contact line and the related problem of 
"true" and "apparent" dynamic contact angles jl]]^] . The 
answer to the riddle must be, in fact, physico-chemical 
rather than purely hydrodynamic, since it depends on 
processes in the immediate vicinity of the three-phase 
boundary. The early detected paradox of a logarithmi- 
cally divergent force required to displace the contact line 
H directly follows from the multivaluedness of the ve- 
locity field at the contact line - if standard viscous hy- 
drodynamics with a no-slip condition on the solid surface 
is to be believed. This paradox has been swept under a 
carpet rather than resolved by introducing a boundary 
condition with a stress or velocity dependent slip 
A drawback of hydrodynamic slip theories lies in their 
inherent inability to predict the dynamic contact angle. 
As a remedy, empirical relationships between the veloc- 
ity and contact angle have to be introduced in model 
computations. 

Clearly, intermolecular forces, that determine the 
static contact angle to begin with, should have a say 
in a dynamic situation. Their direct action is, however, 
restricted to an immediate vicinity of the contact line, 



which is unobservable under available experimental res- 
olution, so that an apparent contact angle seen at meso- 
scopic distances has to be strongly influenced by outer 
hydrodynamic conditions. Near the contact line itself, 
the properties of the fluid are different from those in the 
bulk, and even a common continuum description becomes 
questionable. 

Different approaches to description of the fluid motion 
in the vicinity of the three-phase boundary have been 
tried during the last two decades. The most straight- 
forward way is to introduce intermolecular forces into 
the hydrodynamic equations of motion. This would lead, 
strictly speaking, to very difficult nonlocal equations, in- 
corporating also the effects of variable density and diffuse 
interfaces Even in the sharp interface limit, a non- 
local dependence on the shape of the free interface leads 
to integro-diffcrcntial equations which nobody as yet at- 
tempted to solve. A rational formulation is possible in 
lubrication approximation when the action of inter- 
molecular forces reduces to a simple expression for dis- 
joining pressure between parallel vapor-liquid and liquid- 
solid interphase boundaries Q). This, however, does not 
eliminate the stress singularity, unless in the case of com- 
plete wetting when a sharp contact line is replaced by a 
gradual transition from a precursor film to a liquid film 
of macroscopic thickness 0. At the same time, the usual 
expression for London-van der Waals forces leads to dis- 
joining pressure divergent at small distances and neces- 
sitating a molecular-scale cut-off and leaving the "true" 
contact angle undetermined. This may be formally cor- 
rected by taking account of surface inclination || , but the 
correction becomes effective at non-physical submolecu- 
lar distances. 

A radical solution is abandoning the continuum ap- 
proach altogether in the immediate vicinity of the contact 
line. Slip is feasible on a microscopic scale where it may 
follow from activated diffusion of a first molecular layer 
[^J. Direct numerical simulations of molecular dynamics 
clearly demonstrate the effects of a diffuse boundary and 
effective slip at molecular distances [[l0 11 . Such simula- 
tions, however, cannot involve macroscopic volumes, and 
no ways to incorporate them in a macroscopic descrip- 
tion are known. An alternative approach is to retain 
continuum description but to treat either vapor-liquid, 
or fluid-solid interface, or both as a separate phase with 



1 



properties different from the bulk fluid. This approach 
was adopted by Shikhmurzaev |Q who relied also on de- 
viation from thermodynamic equilibrium near the con- 
tact line as well as on the presence of a residual film to 
avoid the divergences and explain the difference between 
the static and dynamic contact angles. 

Treating the vapor-liquid interface as a separate phase 
may be, indeed, justified when surfactants are present, 
but otherwise a more natural way to account for its spe- 
cial properties of is to consider it as a region interpolating 
between the two phases. The origin of this approach is in 
the diffuse interface model going back to van der Waals 
himself fl"3j |. Much later, it became prominent in the 
phase field models |fl4|| , used mostly in phenomenologi- 
cal theory of solidification where a fictitious phase field, 
rather than density, was used as a continuous variable 
changing across the interphase boundary. The theory of 
van der Waals was widely used for description of equi- 
librium fluid properties, including surface tension and 
line tension in three-phase fluid systems |fl5|| . Applica- 
tions of this theory to dynamical processes in fluids is 
much more difficult, as it requires coupling to hydrody- 
namics. The applicable equations were formulated rather 
recently [ [l6||T^ ]. Seppecher |Q and Jacqmin solved 
the equations of the continuous density field coupled to 
the Stokes equation numerically in a small inner region 
near the contact line, matching it to the outer region 
where the standard sharp-boundary hydrodynamic limit 
applies. The prominent feature of the flow in the inner re- 
gion was a substantial advective mass transport through 
the interphase boundary which served as an effective slip 
mechanism relieving the viscous stress singularity. 

The aim of this communication is a rational analysis 
of the hydrodynamic phase field (diffuse interface) model 
based on the lubrication approximation. After formu- 
lating the basic equations in Section H, we reiterate the 
equilibrium relations defining the surface ten sion o n all 
three kinds of interphase boundaries (Section III A ) and 
discuss appro priate boundary conditions on the solid sur- 
face (Section MB ). This is followed by approximate 
computation of the density profile, equilibrium c hemi- 
cal po tential and energy of the fluid layer (Sections III C , 
HID ). The results of Id computations further serve as 
a basic "vertical" structure of the lubrication theory of 
Section IV, where a slow dependence on the "horizon- 
tal" coordinate is added. We show that the equations 
give a correct sharp-interface limit in both static and dy- 
namic situation. The evolution equation derived in the 
lubrication approximation can be integrated numerically 
yielding the dependence of the spreading velocity on a 
driving force. We shall see in Section IV D that no singu- 
larities develop in the case when the boundary condition 
fixes a unique fluid density at the solid surface. 

This "quasi-equilibrium" theory is modified in Section 
^ where a change of chemical potential across the fluid 
layer is taken into account. We start with discussing 



the "vertical" structure of chemical potential associated 
with viscous and kinetic retardation of steady motion of 
a vapor-liquid interface, and identify the dilute (vapor) 
phase as the locus of substantial variation of chemical po- 
tential. The potential drop is then computed numerically 



in Section VB, yielding a relation between the disjoining 



potential and the flux across isodensity lines. This flux, 
which may be interpreted as incipient evaporation or con- 
densation may help to alleviate viscous stress singularity 
when the boundary conditions make a sharp three-phase 
contact line necessary. 



II. BASIC EQUATIONS 

A general phase field model coupled to hydrodynamics 
includes the following elements: (1) a dynamic equation 
of the phase field variable(s) derived from an appropri- 
ate energy functional; (2) a constituent relation defining 
the dependence of pressure or chemical potential on the 
phase variable(s); (3) the continuity equation; (4) the 
equation for a flow field u(x, t). 

In a one-component system, the appropriate phase 
field variable is density p, General hydrodynamic equa- 
tions for non-equilibrium systems with diffuse interphase 
boundaries are found in the recent review by Anderson 
et al Q. 

The equation for the static density distribution is de- 
rived from the energy functional 



T = J £d 3 x, C = pf( P ) + \K |Vp| ; 



(1) 



where p, is the Lagrange multiplier (chemical potential) 
that serves to insure the mass conservation condition. 
The corresponding Eulcr - Lagrange equation is 



ifV 2 p- W(p))+ M = 0, 



(2) 



We shall suppose that the function f(p) is such that 
Eq. (||) admits two stable solutions, p — p v and p = pi, 
separated by an unstable solution p = p u and p v < p u < 
p%. The solutions are at Maxwell construction (i.e. have 
equal energy pf(p)) at p — 0, so that the chemical po- 
tential can serve as a bias parameter. 

The density field is coupled to hydrodynamics through 
the capillary tensor 



T = Cl-Vp®d£/dVp, 



(3) 



where I is the unity tensor. Eliminating the Lagrange 
multiplier with the help of Eq. (||) yields 

T = {±K\Vp\ 2 + KpV 2 p - p)I - KVp ® Vp, (4) 

where the thermodynamic pressure is defined as p — 
P 2 f(p)- 
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Neglecting the inertial effects, the flow is described by 
the generalized Stokes equation 

V-(T + S) + F = 0, (5) 

where F = — W is an external force and S is the viscous 
stress tensor with the components 

S 3k = riidjVh + d kVj ) + (C - §T7)<5 ife V • v, (6) 

where 77, £ are dynamic viscosities (generally, dependent 
on p), and Vj are components of the velocity field v. The 
system of equations is closed by the continuity equation 

p t + V • (pv) = 0, (7) 

The Stokes equation (g) is rewritten using Eq. (J3J) as 

-\7{p + V) + KpVV 2 p + V - S = 0. (8) 

A more transparent equivalent form, which can be ob- 
tained directly from Eq. (|3|), includes, instead of pressure, 
the chemical potential defined by Eq. (Q): 

- W - KpVp + V • (r/Vv) + V[(C + 3?y)V • v] = 0. 

(9) 

Further on, we shall compute the density and velocity 
field assuming that the characteristic macroscopic length 
L* of the flow field, as well as the scale of density varia- 
tion in the tangential direction far exceed the character- 
istic thickness (if//*) 1 / 2 of the diffuse interface, where 
/* is a characteristic value of f(p). This "thin interface" 
approximation is apt to break down in the vicinity of 
the contact line, unless it is complemented by the "lu- 
brication" approximation, which assumes a small angle 
between the (diffuse) interphase boundary and the solid 
surface. The applicability of this approximation depends 
as well on the boundary conditions at the solid surface. 

III. EQUILIBRIUM RELATIONS 

A. Surface tension and Young— Laplace relation 

Before approaching our main task of the analysis of 
motion in the vicinity of a three-phase boundary, it is 
necessary to clarify relevant properties of the dynamic 
phase field model for the basic case of a diffuse inter- 
face between semi-infinite phases. For a static interface, 
the phase field determines in a usual way the equilibrium 
surface tension jl3],[l5]]. The standard surface tension is 
defined as the energy per unit area of a flat interface sep- 
arating two semi-infinite phases. Static solutions depen- 
dent only on the coordinate z normal to the interface can 
be easily found by solving Eq. . Rescaling the coordi- 
nate by the characteristic width of the diffuse interface, 
and denoting 



g(p) = d p [pf( P )}, (10) 

we have 

p"(z)-g(p) + p = 0. (11) 

The two static solutions are approached at z — » ±00, and 
the boundary is static at p = 0. 

The interfacial energy is computed most easily by using 
as a dependent variable the distortion energy T = \p 2 z - 
Then Eq. (|ll| ) is rewritten as 

T'(p)-g(p)+p = 0, (12) 

Integrating this proves that the distortion energy equals 
the potential energy at any point: 

\p\ = Pf(p) - PP- (13) 
Using this "virial theorem" , we compute 

/oo fPl 
p 2 z dz = / y/2{pf{p) - np)dp. (14) 
-OO J Pv 

Solid-fluid interactions are characterized by an appro- 
priate boundary condition at the solid surface, as elabo- 
rated below. Generally, the density at the solid surface 
will be different in the vapor or liquid phase; we denote 
the respective values as p sv and p s i. Accordingly, the 
"liquid-solid" or "vapor-solid" surface tension 07 or ai is 
computed, respectively, by replacing one of the integra- 
tion limits in Eq. ( p^ ) by p s \ or p sv . 

In the vicinity of a critical point, the appropriate func- 
tion, restricted to small deviations from the critical den- 
sity p c , is a cubic g{p) = p — p c — (p — p c ) 3 - Since our 
aim is a qualitative description of a system far from crit- 
icality involving the vapor phase with negligible density, 
we shall choose a shifted cubic 

g(p)=p(l-2p)(l-p), (15) 

which is at Maxwell construction at p — 0. Then 
f(p) — -?rp(l — p) 2 and the equilibrium surface tensions 
are computed as 

<7 = / p(l - p)dp = ±, 

Jo 

<n = I p(i - p)dp - |(i - Psi?{i + 2pei), 

1 

a v = I p(l-p)dp= y 2 sv (3-2p sv ). (16) 
Jo 

The first formula can be also obtained directly using the 
standard kink solution that approaches p v = at z — > 00 
and pi = 1 at z — > —00: 

p (z) = (l + e*y 1 . (17) 

This solution may be, however, distorted in the vicinity 
of a solid wall. 
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The expressions for a and cr± combine to the Young- 
Laplace formula 



fjj — &i = d cos 8, 



(18) 



where 9 is the "standard" contact angle that should be 
observed at distances much larger than the thickness of 
the transition layer, i.e. unity in the dimensionless units 
ofEq. (0). 

The Young-Laplace formula is a consequence of the 
Noether theorem applied to solutions of Eq. (g) . Suppose 
that the solid surface is coincident with the x axis [z — 0) 
and p(x,z) tends to p g at z — ► oo. Very far on the left 
(x —>■ — oo) the vapor is close to the solid, so that p(x, z) 
tends to a solution p(z) of Eq. ( |il"| ) with p, = 0, such 
that p(0) = p sv and p — ► p v as z — > oo. On the other 
end, for x — ► +oo, the liquid is close to the solid, that 
is p(x,z) tends for x large positive and z <C a; toward 
a solution p(z) of Eq. Jll| ) with = and p(0) = p s i 
and p — > pi as z becomes very large. For a given x, 
there is, however, a value of z, close to tan#, such that 
there is a liquid- vapor interface and for z 3> xtan#, p 
becomes very close to p v , as requested. If p(x, z) satisfies 
these conditions, the liquid-vapor interface is inclined at 
the angle 9 to the solid on scales much larger than the 
microscopic interface thickness (although this angle may 
change at a closer approach). 

The Young-Laplace formula follows from the invari- 
ance of the problem with respect to translations in the 
x direction. Multiply Eq. (|^) (with K rescaled to unity) 
by dp/dx, and integrate over z from z = to oo. This 
yields, after integrating by parts 



d 

dx 



'Mpl-pl)-pf(p)] dz \ = ° 



(19) 



The braced expression is constant along the x axis. This 
constant can be computed for x very large negative and 
very large positive in the configuration just described. 
Equating the results, one gets: 



o; 



;i(^(C)) 2 (cos 2 -sin 2 6) -pf(p)] 



cos 9 ' 



where £ is the coordinate normal to the vapor-liquid in- 
terface, and pd(C) ^ s ^ ne standard kink solution. The 
algebraic term reduces to ^(Po(0) 2 by Eq. (|l3|), and the 



final result is the Young-Laplace formula ( 18 ) . The result 
is not influenced by possible deviations from the standard 
contact angle at a close approach to the solid surface, but, 
of course, hinges on the applicability of Eq. (Q). Since the 
actual inclination angle is apt to change at large distances 
due to external forces, such as gravity or dynamic pres- 
sure, the "standard" angle may be in fact unobservable 
at either small or large distances from the solid. 



B. Boundary conditions 

If the action of the solid on the density field is short- 
range (compared to the thickness of the diffuse interface) , 
it can be accounted for by appropriate boundary condi- 
tions at the solid surface. The boundary conditions are 
usually assigned with the help of the Cahn construction 
|l4],p| balancing the distortion energy, distributed over 
a layer of same order of magnitude as the thickness of 
the diffuse interface, and the energy of fluid-solid inter- 
action concentrated at the boundary. A more consistent 
way to arrive at the same boundary condition is to al- 
low a non-vanishing variation of the density at the solid 
boundary 5p s when the energy functional (Q) is varied. 
In one dimension, this leaves, after integrating by parts, 
the boundary term p'(0)5p s . If the dependence of the 
fluid-solid interaction energy on the fluid density near 
the wall is expressed by a quadratic polynomial 

l(Ps) = 7o - HPs + \lip\-, (20) 

the coefficient at Sp s vanishes, provided 

71-72^ + ^(^)1^=0. (21) 

which is the boundary condition equivalent to that ob- 
tained through the Cahn construction (although the lat- 
ter is expressed in an awkward integral form, including a 
radical with an indefinite sign). 

If one assumes that the solid-fluid interaction is short- 
range compared to the thickness of the diffuse vapor- 
liquid interface, it is likely prevail locally in the vicinity 
of a solid wall. This corresponds to the limiting case 
of very large 71 , 72 , when a simpler Dirichlct boundary 
condition p — p s is enforced on the solid surface. The 
range p v < p s < pi corresponds then to partial wetting. 

With the latter boundary condition and the cubic g(p), 
the contact angle is cos 9 = — l+6p^— 4p 3 ,and is close to 
or 7r when p s is close, respectively, to 1 or —0. If p s = 1— a 
with < a <C 1, we have 9 = 2\/3a. The contact angle 
is zero (complete wetting) at p s > 1. This "standard" 
angle has nothing to do with a "true" contact angle at 
the solid surface. The later is not defined at all in the 
diffuse interface theory, since different isodensity levels 
behave in a qualitatively different way as the solid surface 
is approached. The only level that hits the solid surface 
at the right angle is p = p s ; the levels with p < p s are 
asymtotically parallel, and those with p > p s antiparallcl 
to the surface. 

A more consistent way to derive the boundary condi- 
tion is to start with a general expression for the energy 
of molecular interactions 



T = 



p(x)p(x , )y(|x-x'|)d 3 xd 3 x' 



(22) 



The mean-field energy functional (Q) can be obtained 
from ( p2| ) assuming that the density changes on a char- 
acteristic scale far exceeding the range of the potential 
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V(\k — x'|) and expanding p(x') = p(x) + (x — x') • 
Vp(x) + . . .. The algebraic term in the Lagrangian (|l|) 
is obtained in the zeroth order, and the distortion en- 
ergy in the second order of the expansion. These expres- 
sions are modified when a solid boundary lies within the 
range of the interaction potential. The influence of the 
wall may be particularly strong in the standard case of 
Lennard-Jones interaction potential or a simplified ex- 
pression V oc (x — x')~ 6 with a short-range hard-core 
cut-off, which gives the interaction energy diverging as 
z~ 3 with the distance from the solid surface. The di- 
verging part of the energy may be taken as the surface 
energy potential that has to be minimized to obtain the 
density at the solid surface p s . Under conditions when 
the bulk potential has two minima corresponding to low 
(vapor) and high (liquid) densities, the surface potential 
may also have two minima but the respective values, say, 
p a i and p sv , would be, generally, different from the bulk 
values pi and p v . This brings us to a Dirichlct boundary 
condition similar to that postulated above, but with the 
essential difference that two distinct values are allowed, 
and are likely to be chosen at the solid surface contacting, 
respectively, the liquid and vapor phase. Unlike the case 
when the surface density is unique, all isodensity levels 
in the range p sv < p < p s i hit the solid surface. 

The boundary condition ( pi] ) also allows distinct den- 
sity levels p sv 7^ psi in the areas of the solid surface bor- 
dering either vapor of liquid. Assuming, for example, 
< 71 = a -C 1 , 72 = 0, we have p sv f» a, p s i ~ 1 + fl- 
it appears, however, quite unnatural that the main term 
in the density expansion fixes the density gradient rather 
than the density itself, so that non-monotonic density 
profiles are forbidden in the above example and, on the 
contrary, enforced when 71 is negative. 

C. Density profile in a thin layer 

The interaction between the solid surface and the in- 
terphase boundary can be computed most easily in the 
case when both surfaces are parallel and normal to the 
z axis. The static solution p{z) can be found by solving 
Eq. ( |TT| ) subject to the appropriate boundary conditions 
at the solid wall. Solving the one-dimensional phase field 
equation in the form Eq. ( |l2| ) is elementary; for the cubic 
g(p), the exact solution is expressed in elliptic functions. 
Finding an approximate solution satisfying the bound- 
ary condition p(0) — p s = 1 — a with |a| -C 1 is, however, 
more elucidating. 

We construct the solution by perturbing a standard 
kink solution po(z — h) centered at z = h, e.g. Eq. JT7j ) 
for the cubic g(p). The actual solution is approximated 
to the zero order by the standard kink only when po(~h) 
is sufficiently close to unity; thus, h must satisfy the con- 
dition h > ln(l/a). The density profile is expanded in 
the small parameter o: 



p = po(z — h) + api(z; h) + (23) 

For the time being, we assume p = 0. Then the first- 
order equation is 

p'({z)+g'{ Po ) Pl =Q, (24) 

subject to the boundary condition 

pi(0) = -I + a- x [l - p (-h)] « -1 + V, (25) 

where ip — a~ 1 e~ h < I. 

Due to the exponential decay of interactions, the cor- 
rection to the zero-order solution is actually of a higher 
order of magnitude everywhere except an 0(ln a -1 ) vicin- 
ity of the wall, where po is close to unity. On this interval, 
Eq. (|2^) can be replaced by the equation with constant 
coefficients 

p'{(z)-p 1 (z)=0 > (26) 
The solution decaying at z — > 00 is 

Pl (z) = -e' z (1 - ^) . (27) 
At a > 0, h > ln(2/a), the combined function 

Pa = Po + api = (1 + eJ- h y X - e- z (a - e^) (28) 

reaches a maximum at z = \ ln(ae h — 1) > (Fig. |l|). 
Such a solution describes a liquid layer sandwiched be- 
tween the vapor and the solid. At smaller values of h, the 
maximum disappears, and the solution can be interpreted 
as a pure vapor phase thickening near the solid wall. The 
same solution applies at a < when the density increases 
at the solid surface, whether it is approached from the 
liquid phase or directly from the vapor phase. The ap- 
proximation breaks down at h < ln(l/a), which is, in 
fact, below the minimal possible thickness of the dense 
layer in this model. 

p 




2 4 6 8 10 12 Z 

FIG. 1. Stationary density profiles. Numbers indicate the 
values of the nominal thickness h. 
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If the boundary condition allows two alternative fluid 
densities, a solution with p(0) = a v -C 1 may be also 
possible. This solution, corresponding to vapor phase 
adjacent to the solid surface, is simply p « a v e~ z ; this 
solution can be viewed as a tail of the basic kink centered 
at z = lna„ < (i.e. in the non-physical region). 

Non-monotonic density profiles are unstable. Since, 
however, the influence of the wall decays exponentially 
with the distance, the dynamics is practically frozen 
whenever the interphase boundary is separated from the 
wall by a layer thick compared to the characteristic width 
of the diffuse interface. 



D. Equilibrium chemical potential and energy 

A static solution with a fixed h exists only at a certain 
fixed value of p, which can be determined using a solv- 
ability condition of the first-order equation. In a wider 
context, an appropriate solvability condition serves to ob- 
tain an evolution equation for the nominal position h of 
the interphase boundary. The technique of derivation 
of solvability conditions for a problem involving a semi- 
infinite region and exponentially decaying interactions is 
non-standard and therefore deserves special attention. 

An inhomogeneous first-order equation has a general 
form 



Cpi +H(z) = 0, 



(29) 



containing an inhomogeneity 'H(z) and the linear opera- 
tor 



C = 33 + 9' (Pa 



(30) 



When Eq. (29) is defined on the infinite axis, the solv- 
ability condition of Eq. (j|^) appears due to the presence 
of an eigenfunction of C with zero eigenvalue related to 
the translational symmetry of the kink. The eigenfunc- 
tion, obtained by applying the symmetry operator d/dz, 
is simply p' (z). The solvability condition is fixed by the 
orthogonality of the inhomogeneity to this eigenfunction: 



p' Q {z)H{z)dz = 0. 



(31) 



In the presence of a solid boundary, a difficulty arises, 
however, since the translational invariance is broken and 
no easily computable eigenfunction is available. In addi- 
tion, the orders of magnitudes in the perturbative scheme 
should be estimated in a non-standard way in view of the 
exponential decay of interactions. 

The difficulties are overcome with the help of asymp- 
totic matching technique similar to that employed in the 
theory of vortex dynamics pit] . The solvability condition 



is computed, similar to Eq. (31), using the translational 
eigenfunction on the infinite axis, but the integration is 



not carried out over the entire axis (which now extends 
into the unphysical region z < 0), but starts at some lo- 
cation z = zq > where p differs from the asymptotic 
value p = 1 by an O(a) increment. This generates bound- 
ary terms in the solvability condition, which takes now 
the form 



dpo(z — h) 



dp Q (z 



dz 



dz 

- h) dpi(z) 
dz 



7i(z)dz 



d 2 p (z - h) 
dz 2 



(32) 



z=z 



The boundary values of the first-order solution p\ (z) 
are obtained by solving the first-order equation ( |29|) di- 
rectly on the interval < z < zq, where Eq. (|29| ) can be 
replaced by the equation with constant coefficients (|2q) 
with the added inhomogeneity H(z). The solution of this 
equation is 



where p^ is given by Eq. fl27p and Q{z — () is Green's 
function of Eq. (^6|) . The last term can be neglected for 
certain inhomogeneities, provided the lower limit of the 
integral in the left-hand side of Eq. (|32|) can be shifted to 
— oo without introducing a significant error. The match- 
ing is successful when Eq. (|3^) reduces to a form inde- 
pendent of zq in the leading order. 

The simplest application of the above matching tech- 
nique is the computation of a constant value of chemical 
potential p = p c required to keep the kink at equilib- 
rium (possibly, unstable) at a given location z — h. In 
this case, the inhomogeneity in Eq. (^6| ) is just a constant 
Tt = p c , and the integral in the left-hand side of Eq. ([32] ) 
is p c [po(oo) — po(zo)] = — Pc + 0(a). Since this expres- 
sion remains unchanged in the leading order when zq is 
shifted to — oo, i.e. po(zo) = 1 — 0(a) replaced by unity, 
it is sufficient to use in Eq. © the first term of Eq. (§|) 
only. Retaining the leading term only, we obtain 



Pi (z) 



G(z-QH(QdC, 



(33) 



Pc 



a 2 M(h) = 2a 2 ?A(l - ip) = 1e~ h (a - e~ h ) . (34) 



The first expression demonstrates that the computed 
chemical potentials in fact at most of 0(a 2 ), although 
the equation is nominally of the first order. The gained 
order of magnitude is due to the fast decay of interac- 
tions. Since the computed value is of a higher order, 
there is no need to correct the equilibrium profile com- 
puted in the preceding subsection to 0(a). For a > 0, 
the function p c (h) passes a maximum at the same value 
h = ln(2/a) = O(l) that marks the transition from 
monotonic to non-monotonic density profiles. Sustain- 
ing a static profile requires a bias in favor of the liquid 
state, and the value of p c at the maximum represents the 
critical value of chemical potential required to nucleate a 
thick liquid layer on the solid surface. For a < 0, p c in 



G 



Eq. ( p4| ) is negative and increases monotonically with h; 
in this case, on the contrary, a bias in favor of the vapor 
phase is necessary to keep the interface stationary. 

For the boundary condition p'(0) = —a, the chemical 
potential of the dense solution is equivalent to Eq. ( pi ) 
with the inverted sign of a. The rescaled value is 



M(h) = -2a~ l e~ h (l + 



-1 -h\ 

a e 



(35) 



The correction to energy, defined as 



E(h) 



[pf(p) + ^p 2 z ]dz = a + a 2 E(h), (36) 



also turns out to be of 0(a 2 ). The "virial theorem" used 
in Eq. ( |l4| ) does not hold to this order when density is 
defined by the first-order function p a . The best way to 
compute the energy is to use directly the variational for- 
mulation to relate it with the computed chemical po- 
tential. Requiring the one-dimensional energy functional 
(|l|) to be extremal with respect to h and using in the last 
term p(z) = po(z — h), we compute 



dE 
~dh 



-M I p' (z~h)dz 
lo 



M + 0{a) 



(37) 



IV. MOTION IN A THIN LAYER 

A. Double-scale expansion 

Two-dimensional motion can be rationally treated in 
the familiar "lubrication approximation" , assuming the 
characteristic scale in the "vertical" direction (normal 
to the solid surface) to be much smaller than that in 
the "horizontal" (parallel) direction. When the inter- 
face is weakly inclined and curved, the density is weakly 
dependent on the coordinate x directed along the solid 
surface. Respectively, the vertical velocity v is assumed 
to be much smaller than the horizontal velocity u. The 
scale ratio is determined by the contact angle, and should 
be set at 0(a) = 0(V6) to match the scaling of the 
phase field. The velocities v, u corresponding to weak 
disequilibrium of the phase field considered above will 
be consistently scaled if one assumes d z = 0(1), d x = 
0(VS), u = 0(S 3 ^ 2 ), v = 0(5 2 ). It is further necessary 
for consistent scaling of the hydrodynamic equations that 
the "constant" part of the chemical potential p, associ- 
ated with intcrfacial curvature, disjoining potential, and 
external forces and weakly dependent on x, be of 0(5), 
while the "dynamic" part varying in the vertical direc- 
tion and responsible for motion across isodensity levels, 
be of 0(5 2 ). Further in this Section, we shall assume 
therefore that p + V is independent of z; this assumption 
will be re-examined in Section [v| 



In two dimensions, the term p xx is added to the inho- 
mogeneity in the first-order equation ( p9| ) . In this order, 
the vertical density profile can be represented by the stan- 
dard kink solution po(z — h(x, t)), and the x dependence 
is due to slow variation of h in the "horizontal" direction. 
Thus, 



Pxx = -Poi z - h)h xx + p'o(z - K)h\ 



(38) 



The respective contribution to the solvability condition 
is, in the leading order, 



[p' (z)] 2 dz 



-ah„ 



(39) 



while the contribution of the term containing h 2 vanishes 
in the leading order by symmetry. 

Another possible contribution to the solvability con- 
dition may come from external forces. In the presence 
of gravity directed against the z axis, the equilibrium is 
achieved, according to Eq. (||), at p — po — a 2 Gz rather 
than p = po = const. The rescaled acceleration of grav- 
ity is denoted as a 2 G, which presumes that it matches 
the other terms by the order of magnitude. The integral 
in Eq. ( |j2|) involving the variable part of p is mostly ac- 
cumulated in the diffuse interface region, so that we have 
in the leading order 



G I zp' (z — h)dz « Gh. 
'o 



(40) 



Collecting Eqs. ( |39| ) and (40), we obtain the expression 
for the hydrostatic chemical potential 

p = 6 [M(h) - ah xx + G(h - z)\ , (41) 

where M(h) is defined by Eq. ^ or (||). 



B. Statics in lubrication approximation 

Equation (^) will be used later on to investigate dy- 
namical processes where motion of the contact line is 
involved. In this subsection we investigate the statics of 
this lubrication approximation, and show how it relates 
to the general Young-Laplace result on the static contact 
angle. It seems to be important for the general consis- 
tency of the theory to have dynamical equations for the 
contact angle that reduce to the usual equilibrium the- 
ory in the absence of motion. In most realistic cases, 
the effect of gravity is negligible near the contact line, 
since gravitational forces are much weaker than molecu- 
lar forces. Therefore, the statics of the contact angle, at 
scales in between molecular length scales and the capil- 
lary length (that is the length scale beyond whi ch g ravity 
plays a role), depends on solutions of Equation ( |4l| ) with- 
out the gravity term G(h — z). Moreover, as we want to 
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study equilibrium situations where a liquid-vapor inter- 
face merges with the solid surface, the chemical potential 
p, is set to its equilibrium value, 0, so that the equation 
under consideration is: 



M (h) - ah xx = 0. 



(42) 



In order to derive from this equation the Young- 
Laplace condition, one can use the relation ( |37| ) between 
the energy and chemical potential computed in the end of 
Section HID. We integrate Eq. ( (42| ) with the boundary 
conditions for the function h(x) such that for x — ► — oo, 
the vapor is close to the solid, while for x — > oo, the liq- 
uid is close to the solid, until a height h(x) ~ Ox, 9 <C 1 
where a liquid-vapor interface is situated. The relevant 
first integral of Eq. reads: 



2 ah x 



E{h)-E(h Q ), 



(43) 



where ho is the root of M(h) that gives the thickness of 
the precursor film lying between the solid and the vapor 
phase; ho = ln(l/a) in the model with a cubic f(p) and 
Dirichlet boundary condition. The structure of Eq. fl43| ) 
is obviously similar to the Young-Laplace formula. The 
capillary energy at very large negative x, E{ho), is noth- 
ing but the solid- vapor surface tension, a v . The capillary 
energy at very large positive x where the vapor-liquid in- 
terface is far removed from the solid is the sum of the 
independent contributions of the solid-liquid and a free 
liquid-vapor interfaces. Integrating up to very large pos- 
itive x, where h ~ x9, one gets therefore E(pc) = a + ai. 
Using this in Eq. (^) and subtracting a from both sides 
yields 



a — \oQ 2 rs a cos 8 = a v — 07, 



(44) 



which is the sought after Young-Laplace condition, de- 
rived from the equations of the lubrication approximation 
for the position of the liquid-vapor interface. 



C. Equations of motion in lubrication approximation 

The horizontal velocity u is determined from the hori- 
zontal component of the Stokes equation. Adding gravity 
as an external force, we write the leading order equation 
as 

- po(z - h)P x + (j]u z ) z = 0, (45) 

where the effective pressure P is defined as 

P = Gax + M(h) - ah xx + G(h - z), (46) 

This expression follows from Eq. pi]), with the addition 
of the gravity term acting when the supporting plane 
is weakly inclined. The inclination angle a must be of 



0{\/S) to match by the order of magnitude the other 
terms in the equation. The density profile is given in the 
leading order by the standard kink solution ( |ri| ) centered 
at the nominal interface position h{x) slowly varying in 
the horizontal direction. 

The solution of Eq. (^5|) satisfying the no-slip bound- 
ary condition on the solid boundary and the no stress 
condition at infinity has a general form 



u(z) = rj 1 P x ^(z;h). 



(47) 



The function ^(z;h) depends on an assigned dependence 
of viscosity on density, but the flux upo in the dense layer 
(at z not much larger than h) is nearly the same for ei- 
ther r\ = const or 77 cx p, and is close to the standard 
lubrication solution ^ = —z(h— \z) valid for incom- 
pressible Poiseuillc flow in a layer of thickness h with a 
free boundary. 
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FIG. 2. The function Q(h), compared with the respective 
function for the sharp interface Qo(h) = |/i 3 (dashed line) 

The evolution equation of h is obtained by inserting 
Eqs. (|l7|), (|47]) in the continuity equation (]?]) and inte- 
grating it from to 00. Using the relations 



Ptdz — —ht / p' (z)dz = ht + 0(a), 
Jo 



we obtain 



where 



Q(h) 



(pv) z dz = 0, 



ht - n~ l d x [Q(h)P x ] 



po(z — h)^!{z, h) dz. 



(48) 



(49) 



The function Q(h), computed numerically and plotted in 
Fig. ^, differs only slightly from the respective function 
for the sharp interface Qo(h) — ^h 3 when h exceeds its 
minimal admissible value ho = ln(l/a). Taking into ac- 
count small deviations from the standard kink solution 
near the wall adds only a higher-order correction. 
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D. Quasi-equilibrium spreading 

Apart from a slightly modified volumetric rate, the 
specific contribution of the diffuse interface to Eq. ( p8| ) is 
carried by the function M(h), which is dependent on the 
boundary conditions on the solid surface and expresses 
disjoining potential. It should be emphasized that this 
function is not given a priori but computed in the frame- 
work of the phase field theory (Section [II D| ). The struc- 
ture of Eq. ( fiSj ) is identical to that of standard equations 
of motion of thin liquid films, which are recovered at large 
h when the disjoining potential becomes negligible. At 
small h, the disjoining potential is not singular as in the 
sharp-interface theories with van der Waals interactions 
[PJ. At the same time, the viscous stress singularity at 
the contact line is relaxed as the latter's location becomes 
indefinite. 

Steady flow of a liquid film under the action of disjoin- 
ing potential and gravity can be described by Eq. ( fig| ) 
rewritten in the frame moving with a speed U. We 
shall assume that the liquid layer thickens at x — ► oo, 
and assume U to be positive when the thick layer ad- 
vances. Standard macroscopic arrangements fixing the 
asymptotic conditions at x — > oo are possible, e.g. h — > 
oo, h x = —a for a liquid wedge with the angle a or 
h x = 0, h = y3f7 j olG for an asymptotically flat film on 
an inclined plane. 

Admissible asymptotics at x — > — oo depends on the 
form of the function M(h). If it is given by Eq. (Q) with 
a > 0, the layer may attain asymptotically at x — > -co 
the state of lowest energy h = ho = ln(l/a) (formally, 
this is possible at zero inclination a, although gravity 
effects are negligible in films of molecular thickness) . 

The starting point is Eq. (^) with the effective pres- 
sure given by Eq. ([i"6|). Removing extra parameters by 
rescaling and integrating once yields 



h"'(x) - (M'(h) + G) h'{x) ~aG + 



U(h-h ) 
Q(h) 



0, 



(50) 



where the integration constant has been introduced al- 
lowing for a precursor film with the thickness ho at 
x — > — oo. A more convenient form of Eq. ( |50| ) is ob- 
tained using as the dependent variable y = h 2 , = 2T and 
as the independent variable the nominal thickness h: 



1 



y"(h) - (M'(h) + G) 



J_ f U{h-h ) 
Vv V Q(h) 



aG 



0. 



Equation ( pl| ) is free from singularities which are usu- 
ally caused by divergences of either viscous stress, or dis- 
joining potential, or both, in a layer of vanishing thick- 
ness. It can be integrated numerically starting from the 
asymptotics at x — > -co. The asymptotics of Eq. (|5l|) 



obtained by expanding near h = ho is y >; c 2 (h — ho) 2 , 
implying exponential decay to the "optimal" thickness 
h — ho oc e KX , where the constant k is a positive root of 
the characteristic equation 



k 3 - M'(h )K + U/Q(h ) =0. 



(52) 



Fixing, say, the value of U , one can use the shooting 
method to adjust the value of G satisfying the appropri- 
ate boundary condition at infinity, ^fy = —a. A very 
fine adjustment of the parameter is needed to advance to 
moderate values of h. An example of a computed depen- 
dence of the interface inclination angle on the nominal 
thickness of a dense layer spreading on a horizontal sup- 
port is shown in Fig. ||. 




10 20 30 4C 

FIG. 3. Dependence of the interface inclination angle 9 
on the nominal thickness ft of a spreading dense layer for 
3U — 0.5 and 3(7 = 0.2 (as indicated by numbers at the 
respective curves). The values of G found by shooting are, 
respectively, 0.035123081 and 0.0079817. 



V. NON-EQUILIBRIUM MOTION 

A. Viscously retarded motion 

Equilibrium solutions with p varying along the z axis 
exist only at a particular constant value of /z, equal to 
zero in the adopted gauge. Any deviation of this value 
sets the interface into motion; the interface shift corre- 
sponds to evaporation or condensation retarded by vis- 
cous friction. The simplest case is steady propagation of 
the boundary between two semi-infinite phases. The sta- 
tionary one-dimensional equations in the frame moving 
with the speed c of the steadily propagating interface are 



(pv) z = 0, -pp z + {rjv z ) z = 



(53) 



(51) where v is the single velocity component in this frame; 



external forces are omitted and rj = Q + |tj is the renor- 
malized viscosity, accounting also for the divergence term 
in Eq. (||). These equations are readily integrated yield- 
ing 



j = pv = const, (i = [x c + jR(z), 



(54) 
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where 



R(z) 



ld_ 

p dz 



^dp 



dz 



dz. 



(55) 



The flux j is related to the propagation velocity c as 
j = —c(pi~p v ). The sign of c is chosen in such a way that 
it is positive when the dense (liquid) state advances. The 
constant p c , which may be fixed by external conditions, 
represents the driving force of the process. 

It is reasonable to assume that the disequilibrium is 
weak, so that both p c and the constant flux j are multi- 
plied by a book-keeping small parameter 6 when Eq. ( |5^ ) 
is used in Eq. (|TT|) . The perturbed equation can be ex- 
panded in a usual way, and the relation between the flux 
j and p c is obtained from the solvability condition (pi 



p' (z)R(z)dz. 



(56) 



The integral in the right-hand side can be interpreted 
as the effective friction factor. It depends on the basic 
density profile po(z) as well as on the assumed depen- 
dence of the viscosity on density. If po = p c + P represents 
a weakly perturbed critical density, R{z) = —r] c p~ 3 p z , 
and the integral in Eq. ( p6| ) is proportional to surface 
tension. In the case of vanishing vapor density which 
interests us most, assuming r\ = const leads to a diver- 
gent integral. The divergence is not eliminated also when 
the viscosity is proportional to density. Taking, for ex- 
ample, ff — vp, Eq. (^5|) is evaluated using the relation 
p z = — p(l — p) as R(z) — — i/ha.(po(z)/ p c )- The weak di- 
vergence on the vapor side can be eliminated by assuming 
a small but finite vapor density p v . Then evaluating the 
solvability condition yields 

/oo 
p' (z)R(z)dz 
-CO 

= -cm f' ln^^-dp = -cz/fl + lnpc). (57) 
J Pv Pc 

where p c is the chemical potential at the location with a 
chosen density level p c . 

The dense layer advances (c > 0) at p > 0. This 
causes the chemical potential to drop at at locations with 
lower density ahead of the propagating interface, thereby 
effectively slowing down the advance of the dense layer. 
A sharp drop in the dilute layer, leading to a divergent 
friction factor (p6|), causes substantial deviations from 
the zero-order density profile, which will be taken into 
account in the next section. 



B. Evaporation flux 

We shall consider the case p(oo) — > 0, in view of a 
strong viscous resistance that has lead to the divergence 



of the effective viscosity in Eq. fl5q). Any vertical flux 
causes in this case a substantial change of the chemical 
potential in the vertical direction, as well as a substantial 
distortion of the vertical structure of the density field, as 
will be shown below. 

In a one-dimensional setting, when the flux j is con- 
stant, the vertical structure is computed by solving simul- 
taneously the vertical component of the Stokes equation 
together with Eq. (|Tl|). We assume rj = const (which is 
justified for the dilute phase, where the friction is most 
important) and denote e = fjj. Using p as the indepen- 
dent variable and denoting (p(p) = pi, this equation can 
be rewritten as 



Kp) = g(p) - h<Pp- 



(58) 



This relation can be used in the second Eq. (p3|), yield- 
ing a single equation defining the density profile in the 
presence of evaporation or condensation. The right-hand 
side of this equation can be transformed by replacing 
v z = jdp~ 1 /dz = jp~ 2 y/(p; an apparent change of the 
sign of the last term is due to the fact that p z is negative, 
and has to be defined as —^fip. The resulting equation 
can be integrated once, yielding, after some algebra, 



2 d f f fl( \ 

p T P \Yp- np) 



p 2 



(59) 



The integration constant a can be computed by ap- 
plying this relation deep in the dense layer where ip = p 2 z 
vanishes. This gives a = —pff'(pi) = —p(pi). Thus, 
a is identified with the reverse pressure in the bulk of 
the liquid phase, but is not defined numerically as yet, 
since one has still to compute the shift of the liquid den- 
sity p\ from its standard value pi — 1 due to a shift of 
the chemical potential p. Assuming /i« 1, one can see 
that both pi and a are of the same order of magnitude 
in the bulk of the liquid. One can also observe that, 
as expected, a > at pi < 1 when the dense layer re- 
cedes (evaporates). Deep in the vapor phase, one should 
set p = 0, so that the standard vapor density is not af- 
fected. Nevertheless, if the vapor density tends to zero, 
the boundary condition at the vapor end cannot be ap- 
plied in a straightforward way, since the term containing 
e in Eq. ( |59| ) is indefinite. The solution strategy can be 
outlined then as follows. Picking a certain value of a, we 
integrate Eq. (59) numerically and compute p{pi) using 
Eq. d58|); then a new value of a is computed with the 
help of the algebraic equilibrium relations for the liquid 
phase, and the computation is repeated until it converges 
to a self-consistent solution. 

This procedure can be improved, keeping in mind that 
both p and e are small, though the relation between them, 
crucial for our theory, is still unknown. Equation ( |59| ) 
divided by p 2 can be formally integrated once more and 
rewritten in the form 
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y-pM+p / K{p')d P ' = Q 



where 



K{p) = ep A ^/ip{p) - ap" 



(60) 



(61) 



Differentiating Eq. ( |60| ) yields, in view of Eqs. © and 



p(p) = pK(p) 



K{p')dp'. 



(62) 



The integral accumulates in the "boundary layer" at p — > 
where both terms in Eq. ( |6ll) diverge; thus, computing 
<p(p) in this region is crucial. The scaling in the boundary 
layer is fixed by requiring all terms in Eq. ( |59| ) to be of 
the same order of magnitude. The small parameter e can 
be, indeed, eliminated by setting 



p = re 1 '\ 



(63) 



The rescaled form of Eq. (|60|) applicable in the boundary 
layer is, in the leading order, 



dr \ r 



2\/$ 2 A _ () 



(64) 



Applying the same scaling to Eq. (|62|), one can see, 
however, that the generic estimate is p = 0(e 1 ^ 3 ), which 
is inconsistent with the equilibrium relationships in the 
bulk of the liquid. This can be repaired by adjusting a 
in such a way that the asymptotic value of p at r — > oo 
vanishes in the leading order. This is, indeed, possible, as 
proved by integrating Eq. (Q) numerically. The integra- 
tion starts at some r C using the asymptotic condition 
<f>(r) = A 2 r 4 at r — > 0. A few sample curves p(r) at 
different values of A are drawn in Fig. [IJ all of them ap- 
proach at r — > oo at a certain asymptotic value which has 
to be identified with p(pi). The asymptotic value van- 
ishes atisj 0.677. The residual 0(e 2/3 ) value of p(p t ) 
satisfying the equilibrium relation p(pi) = a, can be ob- 
tained in the next order by allowing an 0(e) deviation of 
a from the chosen value a c ~ 0.677e 2 / 3 . 

Now the solution is completely specified. The chemical 
potential peaks sharply in the dilute phase (Fig. [|). This 
is due to the constraint imposed by a constant flux in one 
dimension, requiring a large driving force in the transi- 
tional layer where a large velocity gradient is necessary 
to compensate the decreasing density. The asymptotics 
dp/dz = —^JTp = — ap 2 corresponds to a rather slow den- 
sity drop-off, p oc z _1 - a dramatic change, compared to 
the exponential decay in Eq. (p"l|). 



-1/3 
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FIG. 4. The chemical potential (i as a function of the 
rescaled density r in the case of evaporation (e > 0). The 
curve corresponding to A — 0.677 with the asymptotic value 
of /i(r) vanishing at r — » oo is flanked by two curves with 
positive and negative values of /x(oo). 

The relation between the chemical potential p(pi) and 
flux j = rj~ 1 e obtained above can be reverted, identifying 
p(pi) = p c with the constant chemical potential in the 
dense layer driving the mass flux j w (0.68?/) "V? /2 - In 
the case of a steadily propagating interface, this can be 
rewritten as a relation between p c and the propagation 
velocity c = — j; unlike Eq. (|56"|), this relation is now 
nonlinear: 



-(0.677frVc /2 - 



C. Condensation flux 



(65) 



The above computation is valid only in the case of 
evaporation (J > or c < 0). A constant condensation 
flux (j < 0) is, clearly, incompatible with vanishing va- 
por density. A positive value of p is required to enable 
condensation, i.e. advance of the dense phase. This, in 
turn, implies a finite vapor density, so that p{p v ) ~ Pv 
at p v <C 1. 

The condensation flux strongly depends on this resid- 
ual density. Equation (|5^) is retained with the sign of 
e inverted, but, since the expression multiplying e is no 
more indefinite at z — > oo, the integration constant a can 
be directly related to p v : 



dp 



2p 



f(p))+plf(Pv) 



0. 



(66) 



The appropriate rescaled variables are again given by 
Eq. (|6^), and the rescaled equation replacing Eq. do4| ) 
reads, in the leading order 



d ( $ \ ( 2 \ <I> 

dr 



^=0 



where (3 
v — ► f3 is 



-1/3 



p v . The asymptotics at 



(67) 



oo or 



11 



$ = n 2 (r- (3f 



4/3 f 



(68) 



This asymptotics corresponds to an exponential decay of 
density to its equilibrium value, r — (3 ~ e KZ at z — + oo. 

Integrating Eq. ( p7[ ) with the asymptotic condition 
(|68|), one can see that also in this case the chemical po- 
tential defined by Eq. ( |58| ) reaches a constant asymptotic 
value at r — > oo . Checking the asymptotics of Eq. ([36]) at 
p —> pi, one can see that the thermodynamic pressure on 
the liquid side, pff'(pi) should be equal that on the vapor 
side, Pvf'(pv), and, hence, be of 0(e 2 / 3 ). This is again 
inconsistent with the generic estimate p — 0(e 1 ' 3 ), im- 
plying, through the equilibrium relationships, the same 
order of magnitude of f'(pi). Hence, as in the preced- 
ing subsection, the value of has to be adjusted in such 
a way that the asymptotic value of p at r — > oo van- 
ishes in the leading order. The value found by shooting 
is P w 0.685 (see Fig. §). 



= 1/3 
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FIG. 5. The chemical potential /i as a function of the 
rescaled density r in the case of condensation (e < 0). The 
curve corresponding to /3 — 0.685 with the asymptotic value 
of p(r) vanishing at r — > oo is flanked by two curves with 
positive and negative values of p(oo). 



D. Kinetically retarded motion 

Taking into account "normal" viscous retardation only 
(with 77 « 77) may exaggerate the actual phase transition 
rate, since transport through a sharp density gradient is 
in fact an activated process, except, perhaps, in an im- 
mediate vicinity of a critical point. When the interface is 
treated as a sharp discontinuity, this may be accounted 
for by introducing a hnite evaporation rate (involving an 
appropriate activation energy) and a condensation "stick- 
ing coefficient" . Both quantities are difficult to estimate 
quantitatively but, in principle, they insure a hnite evap- 
oration or condensation rate even under conditions when 
viscous retardation is absent. 

In the framework of the phase held theory, kinetic re- 
tardation can be accounted for by replacing the station- 
ary equation (Q) or ( |li"| ) by the respective gradient flow 



equation containing a large relaxation time r. In one 
dimension, we have 



rp t = Pzz - g(p) + P- 



(69) 



On the infinite axis, this equation (with p — const) 
has a solution steadily propagating with a speed depen- 
dent on p, and satisfying the stationary equation in the 
comoving frame 



rcp z + p zz - g{p) + p = 0. 



(70) 



In the case of weak disequilibrium, p = 0(5) <C 1, the 
propagation speed c — 0(5) is easily computed, as in the 
preceding subsection, using the solvability condition of 
Eq. ( |70| ) expanded in <5: 



p(pi - p v ) _ 6p 



TG 



(71) 



where a is defined by Eq. (|14| ) and the numerical value 
is given for the cubic g(p)- 

Equations ( |56"| ) and ([7l]) represent two opposite limits 
when, respectively, either viscous or kinetic retardation 
is prevalent. A rough estimate for the lower bound of 
the relaxation time is r oc I 2 /D, where I is the thick- 
ness of the diffuse interface and D is the diffusivity. The 
characteristic time of viscous retardation on the same 
length scale is t v oc l 2 /v, where v = rj/p. For com- 
mon liquids, the Prandtl number Pr = v / D is large, and 
t v /t oc D/v <C 1. Viscous retardation may be still felt 
at larger scales, complementing the kinetic retardation 
near the diffuse boundary. At Pr ^ 1, the flow velocity 
is nearly constant throughout the transitional boundary 
region, and the propagation velocity defined by Eq. ( |7l| ) 
can be viewed as the velocity of the slow drift of the inter- 
phase boundary due to the evaporation or condensation 
in the frame moving with the local velocity of the ambi- 
ent fluid. At fixed propagation velocity, the increments 
due to the viscous and kinetic retardation are additive. 
In the dense layer, the former is negligible at Pr ^ 1, 
although it becomes important in the dilute phase, as we 
have seen in the preceding subsections. 

The scaling of the lubrication approximation (Section 



IV A ) remains consistent only if the relaxation time r in 
Eq. (69) is of 0(8). With this scaling, the speed of the 
vapor-liquid interface displacement is of 0(5 2 ), i.e. of the 
same order of magnitude as the vertical velocity. 



E. Spreading assisted by interphase transport 

The results of the computations that have been carried 
out so far in this section for an infinite fluid layer sepa- 
rated by a diffuse vapor-liquid interface can be applied to 
the spreading problem after minimal modification. In a 
bounded layer, the chemical potential in the liquid phase 
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jj, c driving the evaporation or condensation flux is de- 
termined by the combined action of surface tension and 
disjoining pressure. The disjoining potential can be com- 
puted with the help of the solvability condition, as in Sec- 
tion [II D; the respective formulae remain in force, since 



the flux-related drop of the chemical potential occurs in 
the dilute phase only, and is negligible in the diffuse in- 
terface region, where the translational eigenfunction is 
localized. 

The basic equation of the lubrication approximation, 
Eq. (p8|), is modified in the case of non-equilibrium 
spreading by an added evaporation or condensation term: 



dh 



(72) 



The expressions for the evaporation or condensation flux 
j(P) and even the orders of magnitude vary depending on 
the physical situation under consideration, according to 
the calculations presented in the three preceding subsec- 
tions, with the effective pressure defined by Eq. ( f4(f ) re- 
placing /i c . This determines, in turn, the relative impor- 
tance of the two terms on the right-hand side of Eq. (72). 
In the case of viscously retarded motion with finite p Vl 
j happens to be proportional to ffP, although the term 
representing the horizontal transport through the liquid 
phase is of order rj d 2 P/dx 2 , and is negligible compared 
to the evaporation or condensation term in the lubri- 
cation limit, when the horizontal derivatives are small. 
In this case, flow across the isodensity levels associated 
with evaporation or condensation, driven by the devia- 
tion of the chemical potential from equilibrium, would be 
larger by 0{8~ 1 ) than hydrodynamic "horizontal" mo- 
tion. Therefore, it is likely that this does not represent 
the most usual situations, where evaporation is hindered 
by large activation energies. In the present model a way 
to enter consistently this activation effect is to make the 
evaporation flux and the horizontal transport of the same 
of magnitude. This is done by imposing an 0{8^ 1 ) re- 
laxation time t. Although this connection between a 
molecular quantity r and a macroscopic length scale 6 
may look a bit artificial, this represents a distinguished 
limit where a balance between factors of different physi- 
cal origin is attained. 

The last situation that we have to consider is the one of 
a vapor phase of vanishingly small density when the evap- 
oration flux is related to the jump of chemical potential 
as j oc [ij 2 . In this case, it is possible to have the evapo- 
ration flux and the horizontal transport of the same order 
of magnitude with the choice of scaling r = 0(5^ 1 ^ 2 ). A 
slow density decay caused by evaporation, which might 
lead to a weakly (logarithmically) divergent horizontal 
flux, may be a disturbing factor, but this is certainly an 
artifact caused by a constant flux in a one-dimensional 
setting and not transferable to two-dimensional spread- 
ing. It is of interest to notice at this point that, when j 
is dominant, and in the absence of horizontal flux (which 



would happen far from a solid boundary), one recov- 
ers the classical Thomson expression for the evaporation 
driven by the curvature of the liquid- vapor interface. 

Mass transport across isodensity lines should become 
particularly important when the lubrication approxima- 
tion breaks down. This should happen near the "contact 
line" in the case when two alternative fluid densities near 



the solid wall are possible (see Section [II B). If, say, the 
boundary densities are p sv -C 1 and p s i = 1 — a, a -C 1, 
the three-phase "contact line" can be viewed as a sharp 
transition between O(l) positive and negative values of 
the nominal thickness h, such that e - '' 1 ' <C 1 on cither 
side. This can be treated as a shock of Eq. or (|5Cj). 
The Hugoniot condition which should ensure zero net 
flux through the shock is the equality of chemical poten- 
tials on both sides. Unfortunately, this condition cannot 
be formulated precisely, since the sharp-interface limit 
of the surface tension term is inapplicable in the shock 
region. Moreover, our test computations of the profile 
of the dense layer using Eq. ( |5l| ) with different boundary 
conditions imposed on the "shock" at h = ho showed that 
the spreading velocity is very sensitive to the conditions 
on the shock. 

It remains therefore essential to solve the full system 
of density field and hydrodynamic equation in the shock 
region whenever a sharp transition between alternative 
surface densities is possible. The outer limit of the re- 
solved shock structure should be matched with the lubri- 
cation equations (fl8|), J5C|), or (|72|). The transport across 
isodensity lines in the shock region alleviates the viscous 
stress singularity remaining in the lubrication model. In 
its turn, the latter provides a gradual transition to the 
sharp-interface limit at large distances. 

VI. SUMMARY AND PERSPECTIVES 

As well known, the phase field model provides a sound 
theoretical basis for studying equilibrium capillary phe- 
nomena in fluids. It allows to derive in a straightforward 
manner the classical formulae for the capillary pressure 
and for the equilibrium contact angle, contrary to for- 
mulations based upon the introduction of van der Waals 
forces diverging at short distances. We have shown that 
this model can be extended in a natural way to study a 
thoroughly dynamical spreading process. The lubrication 
limit, where the contact angle is small, allows to derive 
consistently an equation of motion for the liquid-vapor 
interface interacting with the solid surface. In the static 
limit, this equation yields back the equilibrium Young- 
Laplace theory. 

Evaporation or condensation are processes that are in- 
cluded in this model. The driving force for the evapora- 
tion or condensation is the imbalance between the pres- 
sure drop across the interface and its equilibrium value. 
Similarly, and consistently with Seppecher's ps| results, 
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an advancing or receding contact angle differing from its 
equilibrium value makes the contact line a source or sink 
for evaporation or condensation. We suggest to check 
experimentally this interesting phenomenon by observ- 
ing the accumulation of a non-volatile tracer diluted in 
the liquid phase that would be left by evaporation near 
a moving contact line. 

Our analysis indicates that kinetic retardation of the 
interphase transport is essential for a well balanced the- 
ory away from the critical point. The available simula- 
tions of the motion of a diffuse interface near a three- 
phase contact line jl8|,[l9]], taking into account viscous 
retardation only and, in effect, assuming evaporation 
or condensation to be as easy as plain advection, may 
grossly overestimate the rate of interphase transport, but 
the latter remains essential even when its order of mag- 
nitude is reduced due to kinetic retardation. 

The present theory extends itself in a very natural way 
to problems like film breaking. The latter situation is in- 
teresting also because it should allow to approach exper- 
imentally thermodynamical critical points, where phase 
field models certainly apply, although things should be- 
come complicated if the solid-fluid interaction is added 
to the critical phenomena near a moving contact line. 
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